#!/usr/bin/env python
# -*- coding: utf-8 -*-

# AquaCrop crop growth model

import shutil
import subprocess
import datetime
import random
import os
import gc
import re
import math
import sys
import types
import calendar
import glob

import netCDF4 as nc
import numpy as np
import numpy.ma as ma
import pcraster as pcr

import logging
logger = logging.getLogger(__name__)

# file cache to minimize/reduce opening/closing files.  
filecache = dict()

# Global variables:
MV = 1e20
smallNumber = 1E-39

# tuple of netcdf file suffixes (extensions) that can be used:
netcdf_suffixes = ('.nc4','.nc')

# def getFileList(inputDir, filePattern):
# 	'''creates a dictionary of	files meeting the pattern specified'''
# 	fileNameList = glob.glob(os.path.join(inputDir, filePattern))
# 	ll= {}
# 	for fileName in fileNameList:
# 		ll[os.path.split(fileName)[-1]]= fileName
# 	return ll

# def checkVariableInNC(ncFile,varName):

#     logger.debug('Check whether the variable: '+str(varName)+' is defined in the file: '+str(ncFile))
    
#     if ncFile in filecache.keys():
#         f = filecache[ncFile]
#         #~ print "Cached: ", ncFile
#     else:
#         f = nc.Dataset(ncFile)
#         filecache[ncFile] = f
#         #~ print "New: ", ncFile
    
#     varName = str(varName)
    
#     return varName in f.variables.keys()

# def get_dimension_variable(ncFile,dimName):

#     if not checkVariableInNC(ncFile, dimName):            
#         dimvar = None
#     else:
#         if ncFile in filecache.keys():
#             f = filecache[ncFile]
#         else:
#             f = nc.Dataset(ncFile)
#             filecache[ncFile] = f
    
#         dimName = str(dimName)
#         dimvar = f.variables[dimName][:]
    
#     return dimvar

def get_time_index(ncFile, date, useDoy):

    # Get netCDF file and variable name:
    if ncFile in filecache.keys():
        f = filecache[ncFile]
        #~ print "Cached: ", ncFile
    else:
        f = nc.Dataset(ncFile)
        filecache[ncFile] = f
        #~ print "New: ", ncFile
        
    if useDoy == "Yes": 
        logger.debug('Finding the date based on the given climatology doy index (1 to 366, or index 0 to 365)')
        idx = int(dateInput) - 1
    elif useDoy == "month":  
        logger.debug('Finding the date based on the given climatology month index (1 to 12, or index 0 to 11)')
        # make sure that date is in the correct format
        if isinstance(date, str) == True: date = \
                        datetime.datetime.strptime(str(date),'%Y-%m-%d') 
        idx = int(date.month) - 1
    else:
        # make sure that date is in the correct format
        if isinstance(date, str) == True:
            date = datetime.datetime.strptime(str(date),'%Y-%m-%d')
            
        date = datetime.datetime(date.year,date.month,date.day)

        if useDoy == "yearly":
            date  = datetime.datetime(date.year,int(1),int(1))

        if useDoy == "monthly":
            date = datetime.datetime(date.year,date.month,int(1))

        if useDoy == "yearly" or useDoy == "monthly" or useDoy == "daily_seasonal":
            # if the desired year is not available, use the first year or the last year that is available
            first_year_in_nc_file = findFirstYearInNCTime(f.variables['time'])
            last_year_in_nc_file  =  findLastYearInNCTime(f.variables['time'])
            
            if date.year < first_year_in_nc_file:  
                if date.day == 29 and date.month == 2 and calendar.isleap(date.year) and calendar.isleap(first_year_in_nc_file) == False:
                    date = datetime.datetime(first_year_in_nc_file, date.month, 28)
                else:
                    date = datetime.datetime(first_year_in_nc_file, date.month, date.day)
                msg  = "\n"
                msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
                msg += "The date "+str(dateInput)+" is NOT available. "
                msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used."
                msg += "\n"
                logger.warning(msg)
                
            if date.year > last_year_in_nc_file:  
                if date.day == 29 and date.month == 2 and calendar.isleap(date.year) and calendar.isleap(last_year_in_nc_file) == False:
                    date = datetime.datetime(last_year_in_nc_file, date.month, 28)
                else:
                    date = datetime.datetime(last_year_in_nc_file, date.month, date.day)
                msg  = "\n"
                msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
                msg += "The date "+str(dateInput)+" is NOT available. "
                msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used."
                msg += "\n"
                logger.warning(msg)
        try:
            idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, select ='exact')
            msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is available. The 'exact' option is used while selecting netcdf time."
            logger.debug(msg)
        except:
            msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'exact' option CANNOT be used while selecting netcdf time."
            logger.debug(msg)
            try:                                  
                idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, select = 'before')
                msg  = "\n"
                msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
                msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'before' option is used while selecting netcdf time."
                msg += "\n"
            except:
                idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \
                                    select = 'after')
                msg  = "\n"
                msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
                msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'after' option is used while selecting netcdf time."
                msg += "\n"
            logger.warning(msg)

    idx = int(idx)                                                  
    logger.debug('Using the date index '+str(idx))
    return(idx)

def netcdfDim2NumPy(ncFile, dimName, absolutePath = None):
    if absolutePath != None: ncFile = getFullPath(ncFile, absolutePath)
    if ncFile in filecache.keys():
        f = filecache[ncFile]
    else:
        f = nc.Dataset(ncFile)
        filecache[ncFile] = f
    dimName = str(dimName)
    var = f.variables[dimName][:]
    return var
    
def netcdf2PCRobjCloneWithoutTime(ncFile, varName,
                                  cloneMapFileName  = None,\
                                  LatitudeLongitude = True,\
                                  specificFillValue = None,\
                                  absolutePath = None):
    
    if absolutePath != None: ncFile = getFullPath(ncFile, absolutePath)
    
    logger.debug('reading variable: '+str(varName)+' from the file: '+str(ncFile))
    
    # 
    # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file.
    # --- with clone checking
    #     Only works if cells are 'square'.
    #     Only works if cellsizeClone <= cellsizeInput
    # Get netCDF file and variable name:
    if ncFile in filecache.keys():
        f = filecache[ncFile]
        #~ print "Cached: ", ncFile
    else:
        f = nc.Dataset(ncFile)
        filecache[ncFile] = f
        #~ print "New: ", ncFile
    
    #print ncFile
    #f = nc.Dataset(ncFile)  
    varName = str(varName)
    
    if LatitudeLongitude == True:
        try:
            f.variables['lat'] = f.variables['latitude']
            f.variables['lon'] = f.variables['longitude']
        except:
            pass
    
    sameClone = True
    # check whether clone and input maps have the same attributes:
    if cloneMapFileName != None:
        # get the attributes of cloneMap
        attributeClone = getMapAttributesALL(cloneMapFileName)
        cellsizeClone = attributeClone['cellsize']
        rowsClone = attributeClone['rows']
        colsClone = attributeClone['cols']
        xULClone = attributeClone['xUL']
        yULClone = attributeClone['yUL']
        # get the attributes of input (netCDF) 
        cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1]
        cellsizeInput = float(cellsizeInput)
        rowsInput = len(f.variables['lat'])
        colsInput = len(f.variables['lon'])
        xULInput = f.variables['lon'][0]-0.5*cellsizeInput
        yULInput = f.variables['lat'][0]+0.5*cellsizeInput
        # check whether both maps have the same attributes 
        if cellsizeClone != cellsizeInput: sameClone = False
        if rowsClone != rowsInput: sameClone = False
        if colsClone != colsInput: sameClone = False
        if xULClone != xULInput: sameClone = False
        if yULClone != yULInput: sameClone = False
        # print cellsizeClone,cellsizeInput
        # print rowsClone,rowsInput
        # print colsClone,colsInput
        # print xULClone,xULInput
        # print yULClone,yULInput
        # print sameClone

    cropData = f.variables[varName][:,:]       # still original data
    factor = 1                                 # needed in regridData2FinerGrid
    if sameClone == False:
        # crop to cloneMap:
        minX    = min(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)))
        xIdxSta = int(np.where(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)) == minX)[0])
        xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone)))
        minY    = min(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)))
        yIdxSta = int(np.where(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)) == minY)[0])
        yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone)))
        if len(cropData.shape) > 2:
            cropData = f.variables[varName][...,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]
        else:
            cropData = f.variables[varName][yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]
        factor = int(round(float(cellsizeInput)/float(cellsizeClone)))

        if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone)))

    # numpy array
    outnp = regridData2FinerGrid(factor,cropData,MV)

    f = None
    cropData = None 
    return (outnp)
        
def netcdf2NumPyTimeSlice(ncFile,varName,startDate,endDate,
                          useDoy = None,
                          cloneMapFileName = None,
                          LatitudeLongitude = True,
                          specificFillValue = None):

        logger.debug('reading variable: '+str(varName)+' from the file: '+str(ncFile))
    
        if ncFile in filecache.keys():
            f = filecache[ncFile]
            #~ print "Cached: ", ncFile
        else:
            f = nc.Dataset(ncFile)
            filecache[ncFile] = f
            #~ print "New: ", ncFile

        varName = str(varName)

        if LatitudeLongitude == True:
            try:
                f.variables['lat'] = f.variables['latitude']
                f.variables['lon'] = f.variables['longitude']
            except:
                pass

        if varName == "evapotranspiration":        
            try:
                f.variables['evapotranspiration'] = f.variables['referencePotET']
            except:
                pass
        
        # date index
        # NB ncFile should be cached, so don't need to open it again (filecache is a global variable)
        startidx = get_time_index(ncFile, startDate, useDoy)
        endidx   = get_time_index(ncFile, endDate, useDoy)
        idx = np.arange(startidx, endidx + 1)
        
        # compare spatial attributes with clone map
        sameClone = True
        # check whether clone and input maps have the same attributes:
        if cloneMapFileName != None:
            # get the attributes of cloneMap
            attributeClone = getMapAttributesALL(cloneMapFileName)
            cellsizeClone = attributeClone['cellsize']
            rowsClone = attributeClone['rows']
            colsClone = attributeClone['cols']
            xULClone = attributeClone['xUL']
            yULClone = attributeClone['yUL']
            # get the attributes of input (netCDF) 
            cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1]
            cellsizeInput = float(cellsizeInput)
            rowsInput = len(f.variables['lat'])
            colsInput = len(f.variables['lon'])
            xULInput = f.variables['lon'][0]-0.5*cellsizeInput
            yULInput = f.variables['lat'][0]+0.5*cellsizeInput
            # check whether both maps have the same attributes 
            if cellsizeClone != cellsizeInput: sameClone = False
            if rowsClone != rowsInput: sameClone = False
            if colsClone != colsInput: sameClone = False
            if xULClone != xULInput: sameClone = False
            if yULClone != yULInput: sameClone = False

        cropData = f.variables[varName][idx,:,:] # still original data
        factor = 1                               # needed in regridData2FinerGrid

        if sameClone == False:

            logger.debug('Crop to the clone map with lower left corner (x,y): '+str(xULClone)+' , '+str(yULClone))

            # crop to cloneMap:
            #~ xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0])
            minX    = min(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput))) # ; print(minX)
            xIdxSta = int(np.where(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)) == minX)[0])
            xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone)))
            #~ yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0])
            minY    = min(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput))) # ; print(minY)
            yIdxSta = int(np.where(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)) == minY)[0])
            yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone)))
            cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]
            factor = int(round(float(cellsizeInput)/float(cellsizeClone)))
            if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone)))

        # numpy array
        outnp = regridData2FinerGrid(factor,cropData,MV)
        f = None
        cropData = None 
        return (outnp)
                                  
def netcdf2PCRobjClone(ncFile,varName,dateInput,\
                       useDoy = None,
                       cloneMapFileName  = None,\
                       LatitudeLongitude = True,\
                       specificFillValue = None):
    # 
    # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file.
    # --- with clone checking
    #     Only works if cells are 'square'.
    #     Only works if cellsizeClone <= cellsizeInput
    # Get netCDF file and variable name:    
    #~ print ncFile
    
    logger.debug('reading variable: '+str(varName)+' from the file: '+str(ncFile))
    
    if ncFile in filecache.keys():
        f = filecache[ncFile]
        #~ print "Cached: ", ncFile
    else:
        f = nc.Dataset(ncFile)
        filecache[ncFile] = f
        #~ print "New: ", ncFile

    varName = str(varName)

    if LatitudeLongitude == True:
        try:
            f.variables['lat'] = f.variables['latitude']
            f.variables['lon'] = f.variables['longitude']
        except:
            pass
    
    if varName == "evapotranspiration":        
        try:
            f.variables['evapotranspiration'] = f.variables['referencePotET']
        except:
            pass

    # if varName == "kc":   # the variable name in PCR-GLOBWB     
    #    try:
    #        f.variables['kc'] = \
    #             f.variables['Cropcoefficient']  # the variable name in the netcdf file
    #    except:
    #        pass

    # if varName == "interceptCapInput":   # the variable name in PCR-GLOBWB     
    #    try:
    #        f.variables['interceptCapInput'] = \
    #             f.variables['Interceptioncapacity']  # the variable name in the netcdf file
    #    except:
    #        pass

    # if varName == "coverFractionInput":   # the variable name in PCR-GLOBWB     
    #    try:
    #        f.variables['coverFractionInput'] = \
    #             f.variables['Coverfraction']  # the variable name in the netcdf file
    #    except:
    #        pass

    # if varName == "fracVegCover":   # the variable name in PCR-GLOBWB     
    #    try:
    #        f.variables['fracVegCover'] = \
    #             f.variables['vegetation_fraction']  # the variable name in the netcdf file
    #    except:
    #        pass

    # if varName == "minSoilDepthFrac":   # the variable name in PCR-GLOBWB     
    #    try:
    #        f.variables['minSoilDepthFrac'] = \
    #             f.variables['minRootDepthFraction']  # the variable name in the netcdf file
    #    except:
    #        pass

    # if varName == "maxSoilDepthFrac":   # the variable name in PCR-GLOBWB     
    #    try:
    #        f.variables['maxSoilDepthFrac'] = \
    #             f.variables['maxRootDepthFraction']  # the variable name in the netcdf file
    #    except:
    #        pass

    # if varName == "arnoBeta":   # the variable name in PCR-GLOBWB     
    #    try:
    #        f.variables['arnoBeta'] = \
    #             f.variables['arnoSchemeBeta']  # the variable name in the netcdf file
    #    except:
    #        pass

    # date
    date = dateInput
    if useDoy == "Yes": 
        logger.debug('Finding the date based on the given climatology doy index (1 to 366, or index 0 to 365)')
        idx = int(dateInput) - 1
    elif useDoy == "month":  # PS: WE NEED THIS ONE FOR NETCDF FILES that contain only 12 monthly values (e.g. cropCoefficientWaterNC).
        logger.debug('Finding the date based on the given climatology month index (1 to 12, or index 0 to 11)')
        # make sure that date is in the correct format
        if isinstance(date, str) == True: date = \
                        datetime.datetime.strptime(str(date),'%Y-%m-%d') 
        idx = int(date.month) - 1
    else:
        # make sure that date is in the correct format
        if isinstance(date, str) == True: date = \
                        datetime.datetime.strptime(str(date),'%Y-%m-%d') 
        date = datetime.datetime(date.year,date.month,date.day)
        if useDoy == "yearly":
            date  = datetime.datetime(date.year,int(1),int(1))
        if useDoy == "monthly":
            date = datetime.datetime(date.year,date.month,int(1))
        if useDoy == "yearly" or useDoy == "monthly" or useDoy == "daily_seasonal":
            # if the desired year is not available, use the first year or the last year that is available
            first_year_in_nc_file = findFirstYearInNCTime(f.variables['time'])
            last_year_in_nc_file  =  findLastYearInNCTime(f.variables['time'])
            #
            if date.year < first_year_in_nc_file:  
                if date.day == 29 and date.month == 2 and calendar.isleap(date.year) and calendar.isleap(first_year_in_nc_file) == False:
                    date = datetime.datetime(first_year_in_nc_file, date.month, 28)
                else:
                    date = datetime.datetime(first_year_in_nc_file, date.month, date.day)
                msg  = "\n"
                msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
                msg += "The date "+str(dateInput)+" is NOT available. "
                msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used."
                msg += "\n"
                logger.warning(msg)
            if date.year > last_year_in_nc_file:  
                if date.day == 29 and date.month == 2 and calendar.isleap(date.year) and calendar.isleap(last_year_in_nc_file) == False:
                    date = datetime.datetime(last_year_in_nc_file, date.month, 28)
                else:
                    date = datetime.datetime(last_year_in_nc_file, date.month, date.day)
                msg  = "\n"
                msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
                msg += "The date "+str(dateInput)+" is NOT available. "
                msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used."
                msg += "\n"
                logger.warning(msg)
        try:
            idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \
                                select ='exact')
            msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is available. The 'exact' option is used while selecting netcdf time."
            logger.debug(msg)
        except:
            msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'exact' option CANNOT be used while selecting netcdf time."
            logger.debug(msg)
            try:                                  
                idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \
                                    select = 'before')
                msg  = "\n"
                msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
                msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'before' option is used while selecting netcdf time."
                msg += "\n"
            except:
                idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \
                                    select = 'after')
                msg  = "\n"
                msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
                msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'after' option is used while selecting netcdf time."
                msg += "\n"
            logger.warning(msg)
                                                  
    idx = int(idx)                                                  
    logger.debug('Using the date index '+str(idx))

    sameClone = True
    # check whether clone and input maps have the same attributes:
    if cloneMapFileName != None:
        # get the attributes of cloneMap
        attributeClone = getMapAttributesALL(cloneMapFileName)
        cellsizeClone = attributeClone['cellsize']
        rowsClone = attributeClone['rows']
        colsClone = attributeClone['cols']
        xULClone = attributeClone['xUL']
        yULClone = attributeClone['yUL']
        # get the attributes of input (netCDF) 
        cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1]
        cellsizeInput = float(cellsizeInput)
        rowsInput = len(f.variables['lat'])
        colsInput = len(f.variables['lon'])
        xULInput = f.variables['lon'][0]-0.5*cellsizeInput
        yULInput = f.variables['lat'][0]+0.5*cellsizeInput
        # check whether both maps have the same attributes

        # NB: the original code (PCRGLOBWB) used the following line to test
        # equality between cellsize:
        # if cellsizeClone != cellsizeInput: sameClone = False
        # but this is not sensible when dealing with decimal degrees
        # (e.g. 5 arcminute -> 0.08333333) and netCDF files which may have
        # varying levels of precision depending on the creation options of the
        # user. Instead, test almost equality:
        if abs(cellsizeClone - cellsizeInput) > 1e-8: sameClone = False
        if rowsClone != rowsInput: sameClone = False
        if colsClone != colsInput: sameClone = False
        if xULClone != xULInput: sameClone = False
        if yULClone != yULInput: sameClone = False

    cropData = f.variables[varName][int(idx),:,:]       # still original data
    factor = 1                          # needed in regridData2FinerGrid

    if sameClone == False:
        
        logger.debug('Crop to the clone map with lower left corner (x,y): '+str(xULClone)+' , '+str(yULClone))

        # crop to cloneMap:
        #~ xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0])
        minX    = min(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)))
        xIdxSta = int(np.where(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)) == minX)[0])
        xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone)))
        #~ yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0])
        minY    = min(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)))
        yIdxSta = int(np.where(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)) == minY)[0])
        yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone)))
        # cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]
        if len(cropData.shape) > 2:
            cropData = f.variables[varName][idx,...,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]
        else:
            cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]

        factor = int(round(float(cellsizeInput)/float(cellsizeClone)))
        if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone)))

    # numpy array
    outnp = regridData2FinerGrid(factor,cropData,MV)
    
    f = None
    cropData = None 
    return (outnp)

# def netcdf2PCRobjCloneJOYCE(ncFile,varName,dateInput,\
#                        useDoy = None,
#                        cloneMapFileName  = None,\
#                        LatitudeLongitude = True,\
#                        specificFillValue = None):
#     # 
#     # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file.
#     # --- with clone checking
#     #     Only works if cells are 'square'.
#     #     Only works if cellsizeClone <= cellsizeInput
#     # Get netCDF file and variable name:
    
#     #~ print ncFile
    
#     logger.debug('reading variable: '+str(varName)+' from the file: '+str(ncFile))
    
#     if ncFile in filecache.keys():
#         f = filecache[ncFile]
#         #~ print "Cached: ", ncFile
#     else:
#         f = nc.Dataset(ncFile)
#         filecache[ncFile] = f
#         #~ print "New: ", ncFile
    
#     varName = str(varName)
    
#     if LatitudeLongitude == True:
#         try:
#             f.variables['lat'] = f.variables['latitude']
#             f.variables['lon'] = f.variables['longitude']
#         except:
#             pass
    
#     if varName == "evapotranspiration":        
#         try:
#             f.variables['evapotranspiration'] = f.variables['referencePotET']
#         except:
#             pass

#     if varName == "kc":   # the variable name in PCR-GLOBWB     
#        try:
#            f.variables['kc'] = \
#                 f.variables['Cropcoefficient']  # the variable name in the netcdf file
#        except:
#            pass

#     if varName == "interceptCapInput":   # the variable name in PCR-GLOBWB     
#        try:
#            f.variables['interceptCapInput'] = \
#                 f.variables['Interceptioncapacity']  # the variable name in the netcdf file
#        except:
#            pass

#     if varName == "coverFractionInput":   # the variable name in PCR-GLOBWB     
#        try:
#            f.variables['coverFractionInput'] = \
#                 f.variables['Coverfraction']  # the variable name in the netcdf file
#        except:
#            pass

#     if varName == "fracVegCover":   # the variable name in PCR-GLOBWB     
#        try:
#            f.variables['fracVegCover'] = \
#                 f.variables['vegetation_fraction']  # the variable name in the netcdf file
#        except:
#            pass

#     if varName == "arnoBeta":   # the variable name in PCR-GLOBWB     
#        try:
#            f.variables['arnoBeta'] = \
#                 f.variables['arnoSchemeBeta']  # the variable name in the netcdf file
#        except:
#            pass

#     # date
#     date = dateInput
#     if useDoy == "Yes": 
#         logger.debug('Finding the date based on the given climatology doy index (1 to 366, or index 0 to 365)')
#         idx = int(dateInput) - 1
#     else:
#         # make sure that date is in the correct format
#         if isinstance(date, str) == True: date = \
#                         datetime.datetime.strptime(str(date),'%Y-%m-%d') 
#         date = datetime.datetime(date.year,date.month,date.day)
#         if useDoy == "month":
#             logger.debug('Finding the date based on the given climatology month index (1 to 12, or index 0 to 11)')
#             idx = int(date.month) - 1
#         if useDoy == "yearly":
#             date  = datetime.datetime(date.year,int(1),int(1))
#         if useDoy == "monthly":
#             date = datetime.datetime(date.year,date.month,int(1))
#         if useDoy == "yearly" or useDoy == "monthly" or useDoy == "daily_seasonal":
#             # if the desired year is not available, use the first year or the last year that is available
#             first_year_in_nc_file = findFirstYearInNCTime(f.variables['time'])
#             last_year_in_nc_file  =  findLastYearInNCTime(f.variables['time'])
#             #
#             if date.year < first_year_in_nc_file:  
#                 date = datetime.datetime(first_year_in_nc_file,date.month,date.day)
#                 msg  = "\n"
#                 msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
#                 msg += "The date "+str(dateInput)+" is NOT available. "
#                 msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used."
#                 msg += "\n"
#                 logger.warning(msg)
#             if date.year > last_year_in_nc_file:  
#                 date = datetime.datetime(last_year_in_nc_file,date.month,date.day)
#                 msg  = "\n"
#                 msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
#                 msg += "The date "+str(dateInput)+" is NOT available. "
#                 msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used."
#                 msg += "\n"
#                 logger.warning(msg)
#         try:
#             idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \
#                                 select ='exact')
#             msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is available. The 'exact' option is used while selecting netcdf time."
#             logger.debug(msg)
#         except:
#             msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'exact' option CANNOT be used while selecting netcdf time."
#             logger.debug(msg)
#             try:                                  
#                 idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \
#                                     select = 'before')
#                 msg  = "\n"
#                 msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
#                 msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'before' option is used while selecting netcdf time."
#                 msg += "\n"
#             except:
#                 idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \
#                                     select = 'after')
#                 msg  = "\n"
#                 msg += "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!!"+"\n"
#                 msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'after' option is used while selecting netcdf time."
#                 msg += "\n"
#             logger.warning(msg)
                                                  
#     idx = int(idx)                                                  
#     logger.debug('Using the date index '+str(idx))

#     cropData = f.variables[varName][int(idx),:,:].copy()       # still original data
#     factor = 1                                                 # needed in regridData2FinerGrid
    
#     # store latitudes and longitudes to a new variable
#     latitude  = f.variables['lat']
#     longitude = f.variables['lon']
    
#     # check the orientation of the latitude and flip it if necessary
#     we_have_to_flip = False
#     if (latitude[0]- latitude[1]) < 0.0: 
#         we_have_to_flip = True
#         latitude = np.flipud(latitude)
    
#     sameClone = True
#     # check whether clone and input maps have the same attributes:
#     if cloneMapFileName != None:
#         # get the attributes of cloneMap
#         attributeClone = getMapAttributesALL(cloneMapFileName)
#         cellsizeClone = attributeClone['cellsize']
#         rowsClone = attributeClone['rows']
#         colsClone = attributeClone['cols']
#         xULClone = attributeClone['xUL']
#         yULClone = attributeClone['yUL']
#         # get the attributes of input (netCDF) 
#         cellsizeInput = latitude[0]- latitude[1]
#         cellsizeInput = float(cellsizeInput)
#         rowsInput = len(latitude)
#         colsInput = len(longitude)
#         xULInput = longitude[0]-0.5*cellsizeInput
#         yULInput = latitude[0] +0.5*cellsizeInput
#         # check whether both maps have the same attributes 
#         if cellsizeClone != cellsizeInput: sameClone = False
#         if rowsClone != rowsInput: sameClone = False
#         if colsClone != colsInput: sameClone = False
#         if xULClone != xULInput: sameClone = False
#         if yULClone != yULInput: sameClone = False

#     # flip cropData if necessary 
#     if we_have_to_flip: 
#         #~ cropData = cropData[::-1,:]
#         #~ cropData = cropData[::-1,:].copy()

#         #~ cropData = np.flipud(cropData)

#         #~ cropData = np.flipud(cropData)
#         #~ cropData = np.flipud(cropData).copy()

#         #~ original = cropData.copy()
# #~ 
#         #~ print id(cropData)
#         #~ print id(original)

#         #~ cropData = None
#         #~ del cropData
#         #~ cropData = np.flipud(original).copy()
        
#         #~ print type(cropData)
        
#         #~ cropData2 = cropData[::-1,:]
        
#         #~ cropData = None
#         #~ cropData = original[::-1,:]
#         #~ cropData = cropData[::-1,:]

#         cropData = cropData[::-1,:]
        
#         print type(cropData)

#         print "Test test tet"
#         print id(cropData)
#         #~ print id(original)

#         #~ cropData = cropData[::-1,:].copy()

#         pcr_map = pcr.numpy2pcr(pcr.Scalar, cropData, -999.9)
#         pcr.report(pcr_map, "test2.map")
#         os.system("aguila test2.map")
    
#     if sameClone == False:
        
#         logger.debug('Crop to the clone map with lower left corner (x,y): '+str(xULClone)+' , '+str(yULClone))
#         # crop to cloneMap:
#         minX    = min(abs(longitude[:] - (xULClone + 0.5*cellsizeInput))) # ; print(minX)
#         xIdxSta = int(np.where(abs(longitude[:] - (xULClone + 0.5*cellsizeInput)) == minX)[0])
#         xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone)))
#         minY    = min(abs(latitude[:] - (yULClone - 0.5*cellsizeInput))) # ; print(minY)
#         yIdxSta = int(np.where(abs(latitude[:] - (yULClone - 0.5*cellsizeInput)) == minY)[0])
#         yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone)))
#         cropData = cropData[yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]

#         factor = int(round(float(cellsizeInput)/float(cellsizeClone)))
#         if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone)))

#     # convert to PCR object and close f
#     if specificFillValue != None:
#         outPCR = pcr.numpy2pcr(pcr.Scalar, \
#                   regridData2FinerGrid(factor,cropData,MV), \
#                   float(specificFillValue))
#     else:
#         outPCR = pcr.numpy2pcr(pcr.Scalar, \
#                   regridData2FinerGrid(factor,cropData,MV), \
#                   float(f.variables[varName]._FillValue))
                  
#     #f.close();
#     f = None ; cropData = None 
#     # PCRaster object
#     return (outPCR)


# def netcdf2PCRobjCloneWindDist(ncFile,varName,dateInput,useDoy = None,
#                        cloneMapFileName=None):
#     # EHS (02 SEP 2013): This is a special function made by Niko Wanders (for his DA framework).
#     # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file.
#     # --- with clone checking
#     #     Only works if cells are 'square'.
#     #     Only works if cellsizeClone <= cellsizeInput
    
#     # Get netCDF file and variable name:
#     f = nc.Dataset(ncFile)
#     varName = str(varName)

#     # date
#     date = dateInput
#     if useDoy == "Yes": 
#         idx = dateInput - 1
#     else:
#         if isinstance(date, str) == True: date = \
#                         datetime.datetime.strptime(str(date),'%Y-%m-%d') 
#         date = datetime.datetime(date.year,date.month,date.day)
#         # time index (in the netCDF file)
#         nctime = f.variables['time']  # A netCDF time variable object.
#         idx = nc.date2index(date, nctime, calendar=nctime.calendar, \
#                                                      select='exact')
#     idx = int(idx)                                                  

#     sameClone = True
#     # check whether clone and input maps have the same attributes:
#     if cloneMapFileName != None:
#         # get the attributes of cloneMap
#         attributeClone = getMapAttributesALL(cloneMapFileName)
#         cellsizeClone = attributeClone['cellsize']
#         rowsClone = attributeClone['rows']
#         colsClone = attributeClone['cols']
#         xULClone = attributeClone['xUL']
#         yULClone = attributeClone['yUL']
#         # get the attributes of input (netCDF) 
#         cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1]
#         cellsizeInput = float(cellsizeInput)
#         rowsInput = len(f.variables['lat'])
#         colsInput = len(f.variables['lon'])
#         xULInput = f.variables['lon'][0]-0.5*cellsizeInput
#         yULInput = f.variables['lat'][0]+0.5*cellsizeInput
#         # check whether both maps have the same attributes 
#         if cellsizeClone != cellsizeInput: sameClone = False
#         if rowsClone != rowsInput: sameClone = False
#         if colsClone != colsInput: sameClone = False
#         if xULClone != xULInput: sameClone = False
#         if yULClone != yULInput: sameClone = False

#     cropData = f.variables[varName][int(idx),:,:]       # still original data
#     factor = 1                          # needed in regridData2FinerGrid
#     if sameClone == False:
#         # crop to cloneMap:
#         xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0])
#         xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone)))
#         yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0])
#         yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone)))
#         cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]
#         factor = int(float(cellsizeInput)/float(cellsizeClone))
    
#     # convert to PCR object and close f
#     outPCR = pcr.numpy2pcr(pcr.Scalar, \
#                regridData2FinerGrid(factor,cropData,MV), \
#                   float(0.0))
#     f.close();
#     f = None ; cropData = None 
#     # PCRaster object
#     return (outPCR)    
    
# def netcdf2PCRobjCloneWind(ncFile,varName,dateInput,useDoy = None,
#                        cloneMapFileName=None):
#     # EHS (02 SEP 2013): This is a special function made by Niko Wanders (for his DA framework).
#     # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file.
#     # --- with clone checking
#     #     Only works if cells are 'square'.
#     #     Only works if cellsizeClone <= cellsizeInput
    
#     # Get netCDF file and variable name:
#     f = nc.Dataset(ncFile)
#     varName = str(varName)

#     # date
#     date = dateInput
#     if useDoy == "Yes": 
#         idx = dateInput - 1
#     else:
#         if isinstance(date, str) == True: date = \
#                         datetime.datetime.strptime(str(date),'%Y-%m-%d') 
#         date = datetime.datetime(date.year,date.month,date.day, 0, 0)
#         # time index (in the netCDF file)
#         nctime = f.variables['time']  # A netCDF time variable object.
#         idx = nc.date2index(date, nctime, select="exact")
#     idx = int(idx)                                                  

#     sameClone = True
#     # check whether clone and input maps have the same attributes:
#     if cloneMapFileName != None:
#         # get the attributes of cloneMap
#         attributeClone = getMapAttributesALL(cloneMapFileName)
#         cellsizeClone = attributeClone['cellsize']
#         rowsClone = attributeClone['rows']
#         colsClone = attributeClone['cols']
#         xULClone = attributeClone['xUL']
#         yULClone = attributeClone['yUL']
#         # get the attributes of input (netCDF) 
#         cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1]
#         cellsizeInput = float(cellsizeInput)
#         rowsInput = len(f.variables['lat'])
#         colsInput = len(f.variables['lon'])
#         xULInput = f.variables['lon'][0]-0.5*cellsizeInput
#         yULInput = f.variables['lat'][0]+0.5*cellsizeInput
#         # check whether both maps have the same attributes 
#         if cellsizeClone != cellsizeInput: sameClone = False
#         if rowsClone != rowsInput: sameClone = False
#         if colsClone != colsInput: sameClone = False
#         if xULClone != xULInput: sameClone = False
#         if yULClone != yULInput: sameClone = False

#     cropData = f.variables[varName][int(idx),:,:]       # still original data
#     factor = 1                          # needed in regridData2FinerGrid
#     if sameClone == False:
#         # crop to cloneMap:
#         xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0])
#         xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone)))
#         yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0])
#         yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone)))
#         cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]
#         factor = int(float(cellsizeInput)/float(cellsizeClone))
    
#     # convert to PCR object and close f
#     outPCR = pcr.numpy2pcr(pcr.Scalar, \
#                regridData2FinerGrid(factor,cropData,MV), \
#                   float(f.variables[varName]._FillValue))
#     f.close();
#     f = None ; cropData = None 
#     # PCRaster object
#     return (outPCR)    
    
# def netcdf2PCRobj(ncFile,varName,dateInput):
#     # EHS (04 APR 2013): To convert netCDF (tss) file to PCR file.
#     # The cloneMap is globally defined (outside this method).
    
#     # Get netCDF file and variable name:
#     f = nc.Dataset(ncFile)
#     varName = str(varName)

#     # date
#     date = dateInput
#     if isinstance(date, str) == True: date = \
#                     datetime.datetime.strptime(str(date),'%Y-%m-%d') 
#     date = datetime.datetime(date.year,date.month,date.day)
    
#     # time index (in the netCDF file)
#     nctime = f.variables['time']  # A netCDF time variable object.
#     idx = nc.date2index(date, nctime, calendar=nctime.calendar, \
#                                                  select='exact') 
    
#     # convert to PCR object and close f
#     outPCR = pcr.numpy2pcr(pcr.Scalar,(f.variables[varName][idx].data), \
#                              float(f.variables[varName]._FillValue))
#     f.close(); f = None ; del f
#     # PCRaster object
#     return (outPCR)

# def makeDir(directoryName):
#     try:
#         os.makedirs(directoryName)
#     except OSError:
#         pass

# def writePCRmapToDir(v,outFileName,outDir):
#     # v: inputMapFileName or floating values
#     # cloneMapFileName: If the inputMap and cloneMap have different clones,
#     #                   resampling will be done. Then,   
#     fullFileName = getFullPath(outFileName,outDir)
#     logger.debug('Writing a pcraster map to : '+str(fullFileName))
#     pcr.report(v,fullFileName)

def readPCRmapClone(v,cloneMapFileName,tmpDir,absolutePath=None,isLddMap=False,cover=None,isNomMap=False):
	# v: inputMapFileName or floating values
	# cloneMapFileName: If the inputMap and cloneMap have different clones,
	#                   resampling will be done.   
    logger.debug('read file/values: '+str(v))
    if v == "None":
        #~ PCRmap = str("None")
        PCRmap = None                                                   # 29 July: I made an experiment by changing the type of this object. 
    elif not re.match(r"[0-9.-]*$",v):
        if absolutePath != None: v = getFullPath(v,absolutePath)
        # print v
        # print cloneMapFileName
        sameClone = isSameClone(v,cloneMapFileName)
        if sameClone == True:
            PCRmap = pcr.readmap(v)
        else:
            # resample using GDAL:
            output = tmpDir+'temp.map'
            warp = gdalwarpPCR(v,output,cloneMapFileName,tmpDir,isLddMap,isNomMap)
            # read from temporary file and delete the temporary file:
            PCRmap = pcr.readmap(output)
            if isLddMap == True: PCRmap = pcr.ifthen(pcr.scalar(PCRmap) < 10., PCRmap)
            if isLddMap == True: PCRmap = pcr.ldd(PCRmap)
            if isNomMap == True: PCRmap = pcr.ifthen(pcr.scalar(PCRmap) >  0., PCRmap)
            if isNomMap == True: PCRmap = pcr.nominal(PCRmap)
            if os.path.isdir(tmpDir):
                shutil.rmtree(tmpDir)
            os.makedirs(tmpDir)
    else:
        PCRmap = pcr.spatial(pcr.scalar(float(v)))
    if cover != None:
        PCRmap = pcr.cover(PCRmap, cover)
    co = None; cOut = None; err = None; warp = None
    del co; del cOut; del err; del warp
    stdout = None; del stdout
    stderr = None; del stderr

    # SM: revisit this
    PCRmap = pcr.pcr2numpy(PCRmap, np.nan)
    
    return PCRmap    

# def readPCRmap(v):
# 	# v : fileName or floating values
#     if not re.match(r"[0-9.-]*$", v):
#         PCRmap = pcr.readmap(v)
#     else:
#         PCRmap = pcr.scalar(float(v))
#     return PCRmap    

def isSameClone(inputMapFileName,cloneMapFileName):    
    # reading inputMap:
    attributeInput = getMapAttributesALL(inputMapFileName)
    cellsizeInput = attributeInput['cellsize']
    rowsInput = attributeInput['rows']
    colsInput = attributeInput['cols']
    xULInput = attributeInput['xUL']
    yULInput = attributeInput['yUL']
    # reading cloneMap:
    attributeClone = getMapAttributesALL(cloneMapFileName)
    cellsizeClone = attributeClone['cellsize']
    rowsClone = attributeClone['rows']
    colsClone = attributeClone['cols']
    xULClone = attributeClone['xUL']
    yULClone = attributeClone['yUL']
    # check whether both maps have the same attributes? 
    sameClone = True
    if cellsizeClone != cellsizeInput: sameClone = False
    if rowsClone != rowsInput: sameClone = False
    if colsClone != colsInput: sameClone = False
    if xULClone != xULInput: sameClone = False
    if yULClone != yULInput: sameClone = False
    return sameClone

# def gdalwarpPCR(input,output,cloneOut,tmpDir,isLddMap=False,isNominalMap=False):
#     # 19 Mar 2013 created by Edwin H. Sutanudjaja
#     # all input maps must be in PCRaster maps
#     # 
#     # remove temporary files:
#     co = 'rm '+str(tmpDir)+'*.*'
#     cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate()
#     # 
#     # converting files to tif:
#     co = 'gdal_translate -ot Float64 '+str(input)+' '+str(tmpDir)+'tmp_inp.tif'
#     if isLddMap == True: co = 'gdal_translate -ot Int32 '+str(input)+' '+str(tmpDir)+'tmp_inp.tif'
#     if isNominalMap == True: co = 'gdal_translate -ot Int32 '+str(input)+' '+str(tmpDir)+'tmp_inp.tif'
#     cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate()
#     # 
#     # get the attributes of PCRaster map:
#     cloneAtt = getMapAttributesALL(cloneOut)
#     xmin = cloneAtt['xUL']
#     ymin = cloneAtt['yUL'] - cloneAtt['rows']*cloneAtt['cellsize']
#     xmax = cloneAtt['xUL'] + cloneAtt['cols']*cloneAtt['cellsize']
#     ymax = cloneAtt['yUL'] 
#     xres = cloneAtt['cellsize']
#     yres = cloneAtt['cellsize']
#     te = '-te '+str(xmin)+' '+str(ymin)+' '+str(xmax)+' '+str(ymax)+' '
#     tr = '-tr '+str(xres)+' '+str(yres)+' '
#     co = 'gdalwarp '+te+tr+ \
#          ' -srcnodata -3.4028234663852886e+38 -dstnodata mv '+ \
#            str(tmpDir)+'tmp_inp.tif '+ \
#            str(tmpDir)+'tmp_out.tif'
#     cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate()
#     # 
#     co = 'gdal_translate -of PCRaster '+ \
#               str(tmpDir)+'tmp_out.tif '+str(output)
#     cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate()
#     # 
#     co = 'mapattr -c '+str(cloneOut)+' '+str(output)
#     cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate()
#     # 
#     #~ co = 'aguila '+str(output)
#     #~ print(co)
#     #~ cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate()
#     # 
#     co = 'rm '+str(tmpDir)+'tmp*.*'
#     cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate()
#     co = None; cOut = None; err = None
#     del co; del cOut; del err
#     stdout = None; del stdout
#     stderr = None; del stderr
#     n = gc.collect() ; del gc.garbage[:] ; n = None ; del n

def getFullPath(inputPath,absolutePath,completeFileName = True):
    # 19 Mar 2013 created by Edwin H. Sutanudjaja
    # Function: to get the full absolute path of a folder or a file
          
    # replace all \ with /
    inputPath = str(inputPath).replace("\\", "/")
    absolutePath = str(absolutePath).replace("\\", "/")
    
    # tuple of suffixes (extensions) that can be used:
    suffix = ('/','_','.nc4','.map','.nc','.dat','.txt','.asc','.ldd','.tbl',\
              '.001','.002','.003','.004','.005','.006',\
              '.007','.008','.009','.010','.011','.012')
    
    if inputPath.startswith('/') or str(inputPath)[1] == ":":
        fullPath = str(inputPath)
    else:
        if absolutePath.endswith('/'): 
            absolutePath = str(absolutePath)
        else:
            absolutePath = str(absolutePath)+'/'
        fullPath = str(absolutePath)+str(inputPath)
    
    if completeFileName:
        if fullPath.endswith(suffix): 
            fullPath = str(fullPath)
        else:
            fullPath = str(fullPath)+'/'    

    return fullPath    		

# def findISIFileName(year,model,rcp,prefix,var):
#     histYears = [1951,1961,1971,1981,1991,2001]
#     sYears = [2011,2021,2031,2041,2051,2061,2071,2081,2091]
#     rcpStr = rcp
#     if year >= sYears[0]:
#         sYear = [i for i in range(len(sYears)) if year >= sYears[i]]
#         sY  = sYears[sYear[-1]]
        
#     elif year < histYears[-1]:
       
#         sYear = [i for i in range(len(histYears)) if year >= histYears[i] ]
#         sY  = histYears[sYear[-1]]
    
#     if year >= histYears[-1] and year < sYears[0]:
         
#         if model == 'HadGEM2-ES':
#             if year < 2005:
#                 rcpStr = 'historical'               
#                 sY = 2001
#                 eY = 2004
#             else:
#                 rcpStr = rcp
#                 sY = 2005
#                 eY = 2010
#         if model == 'IPSL-CM5A-LR' or model == 'GFDL-ESM2M':
#             if year < 2006:
#                 rcpStr = 'historical'
#                 sY = 2001
#                 eY = 2005
#             else:
#                 rcpStr = rcp
#                 sY = 2006
#                 eY = 2010
            
#     else:        
#         eY = sY + 9
#         if sY == 2091:
#             eY  = 2099
#     if model == 'HadGEM2-ES':
#         if year < 2005:
#             rcpStr = 'historical'
#     if model == 'IPSL-CM5A-LR' or model == 'GFDL-ESM2M':
#         if year < 2006:
#             rcpStr = 'historical'
#     #print year,sY,eY
#     return "%s_%s_%s_%s_%i-%i.nc" %(var,prefix,model.lower(),rcpStr,sY,eY)
 
# def get_random_word(wordLen):
#     word = ''
#     for i in range(wordLen):
#         word += random.choice('ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789')
#     return word
    
# def isLastDayOfMonth(date):
#     if (date + datetime.timedelta(days=1 )).day == 1:
#         return True
#     else:
#         return False

def getMapAttributesALL(cloneMap,arcDegree=True):
    cOut,err = subprocess.Popen(str('mapattr -p %s ' %(cloneMap)), stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate()

    if err !=None or cOut == []:
        print("Something wrong with mattattr in VirtualOS, maybe clone Map does not exist ?")
        sys.exit()
    cellsize = float(cOut.split()[7])
    if arcDegree == True: cellsize = round(cellsize * 360000.)/360000.
    mapAttr = {'cellsize': float(cellsize)        ,\
               'rows'    : float(cOut.split()[3]) ,\
               'cols'    : float(cOut.split()[5]) ,\
               'xUL'     : float(cOut.split()[17]),\
               'yUL'     : float(cOut.split()[19])}
    co = None; cOut = None; err = None
    del co; del cOut; del err
    n = gc.collect() ; del gc.garbage[:] ; n = None ; del n
    return mapAttr 

# def getMapAttributes(cloneMap,attribute,arcDegree=True):
#     cOut,err = subprocess.Popen(str('mapattr -p %s ' %(cloneMap)), stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate()
#     #print cOut
#     if err !=None or cOut == []:
#         print "Something wrong with mattattr in VirtualOS, maybe clone Map does not exist ? "
#         sys.exit()
#     #print cOut.split()
#     co = None; err = None
#     del co; del err
#     n = gc.collect() ; del gc.garbage[:] ; n = None ; del n
#     if attribute == 'cellsize':
#         cellsize = float(cOut.split()[7])
#         if arcDegree == True: cellsize = round(cellsize * 360000.)/360000.
#         return cellsize  
#     if attribute == 'rows':
#         return int(cOut.split()[3])
#         #return float(cOut.split()[3])
#     if attribute == 'cols':
#         return int(cOut.split()[5])
#         #return float(cOut.split()[5])
#     if attribute == 'xUL':
#         return float(cOut.split()[17])
#     if attribute == 'yUL':
#         return float(cOut.split()[19])
    
# def getMapTotal(mapFile):
#     ''' outputs the sum of all values in a map file '''

#     total, valid = pcr.cellvalue(pcr.maptotal(mapFile),1)
#     return total

# def getMapTotalHighPrecisionButOnlyForPositiveValues_NEEDMORETEST(mapFile):
#     ''' outputs the sum of all values in a map file '''

#     # STILL UNDER DEVELOPMENT - NOT FULLY TESTED
    
#     # input map - note that all values must be positive
#     remainingMapValue = pcr.max(0.0, mapFile)
    
#     # loop from biggest values
#     min_power_number = 0                                                      # The minimum value is zero.
#     max_power_number = int(pcr.mapmaximum(pcr.log10(remainingMapValue))) + 1
#     step = 1
#     total_map_for_every_power_number = {}
#     for power_number in range(max_power_number, min_power_number - step, -step):
        
#         # cell value in this loop        
#         currentCellValue = pcr.rounddown(remainingMapValue * pcr.scalar(10.**(power_number))) / pcr.scalar(10.**(power_number))
#         if power_number == min_power_number: currentCellValue = remainingMapValue
        
#         # map total in this loop
#         total_in_this_loop, valid = pcr.cellvalue(pcr.maptotal(currentCellValue), 1)
#         total_map_for_every_power_number[str(power_number)] = total_in_this_loop
                
#         # remaining map value 
#         remainingMapValue = pcr.max(0.0, remainingMapValue - currentCellValue)
        
#     # sum from the smallest values (minimizing numerical errors)
#     total = pcr.scalar(0.0)
#     for power_number in range(min_power_number, max_power_number + step, step):
#         total += total_map_for_every_power_number[str(power_number)]

#     return total

# def get_rowColAboveThreshold(map, threshold):
#     npMap = pcr.pcr2numpy(map, -9999)
#     (nr, nc) = np.shape(npMap)
#     for r in range(0, nr):
#         for c in range(0, nc):
#             if npMap[r, c] != -9999:
#                 if np.abs(npMap[r, c]) > threshold:


#                     return (r, c)

def getLastDayOfMonth(date):
    ''' returns the last day of the month for a given date '''

    if date.month == 12:
        return date.replace(day=31)
    return date.replace(month=date.month + 1, day=1) - datetime.timedelta(days=1)



# def getMinMaxMean(mapFile,ignoreEmptyMap=False):
#     mn = pcr.cellvalue(pcr.mapminimum(mapFile),1)[0]
#     mx = pcr.cellvalue(pcr.mapmaximum(mapFile),1)[0]
#     nrValues = pcr.cellvalue(pcr.maptotal(pcr.scalar(pcr.defined(mapFile))), 1 )[0] #/ getNumNonMissingValues(mapFile)
#     if nrValues == 0.0 and ignoreEmptyMap: 
#         return 0.0,0.0,0.0
#     else:
#         return mn,mx,(getMapTotal(mapFile) / nrValues)

# def getMapVolume(mapFile, cellareaFile):
#     ''' returns the sum of all grid cell values '''
#     volume = mapFile * cellareaFile
#     return (getMapTotal(volume) / 1)

# def secondsPerDay():
#     return float(3600 * 24)
    
# def getValDivZero(x,y,y_lim=smallNumber,z_def= 0.):
#   #-returns the result of a division that possibly involves a zero
#   # denominator; in which case, a default value is substituted:
#   # x/y= z in case y > y_lim,
#   # x/y= z_def in case y <= y_lim, where y_lim -> 0.
#   # z_def is set to zero if not otherwise specified
#   return pcr.ifthenelse(y > y_lim,x/pcr.max(y_lim,y),z_def)

# def getValFloatDivZero(x,y,y_lim,z_def= 0.):
#   #-returns the result of a division that possibly involves a zero
#   # denominator; in which case, a default value is substituted:
#   # x/y= z in case y > y_lim,
#   # x/y= z_def in case y <= y_lim, where y_lim -> 0.
#   # z_def is set to zero if not otherwise specified
#   if y > y_lim:
#     return x / max(y_lim,y)
#   else:
#     return z_def

# def retrieveMapValue(pcrX,coordinates):
#     #-retrieves values from a map and returns an array conform the IDs stored in properties
#     nrRows= coordinates.shape[0]
#     x= np.ones((nrRows))* MV
#     tmpIDArray= pcr.pcr2numpy(pcrX,MV)
#     for iCnt in xrange(nrRows):
#       row,col= coordinates[iCnt,:]
#       if row != MV and col != MV:
#         x[iCnt]= tmpIDArray[row,col]
#     return x

# def returnMapValue(pcrX,x,coord):
#     #-retrieves value from an array and update values in the map
#     if x.ndim == 1:
#       nrRows= 1

#     tempIDArray= pcr.pcr2numpy(pcrX,MV)
#     #print tempIDArray
#     temporary= tempIDArray
#     nrRows= coord.shape[0]
#     for iCnt in xrange(nrRows):
#       row,col= coord[iCnt,:]
#       if row != MV and col != MV:
#         tempIDArray[row,col]= (x[iCnt])
#        # print iCnt,row,col,x[iCnt]
#     pcrX= pcr.numpy2pcr(pcr.Scalar,tempIDArray,MV)
#     return pcrX
    
# def getQAtBasinMouths(discharge, basinMouth):
#     temp = pcr.ifthenelse(basinMouth != 0 , discharge * secondsPerDay(),0.)
#     pcr.report(temp,"temp.map")
#     return (getMapTotal(temp)  / 1e9)

# def regridMapFile2FinerGrid (rescaleFac,coarse):
#     if rescaleFac ==1:
#         return coarse
#     return pcr.numpy2pcr(pcr.Scalar, regridData2FinerGrid(rescaleFac,pcr.pcr2numpy(coarse,MV),MV),MV)

def regridData2FinerGrid(rescaleFac, coarse, MV):
    # TODO: check with >3 dimensions (only checked with 2D
    # and 3D arrays so far)
    
    if rescaleFac == 1:
        return coarse

    dim = np.shape(coarse)
    ndim = len(dim)
    nr,nc = dim[-2:]
    if ndim > 2:
        othdim = dim[0:(ndim-2)]
        othdimprod = reduce(mul, othdim)
        fine = np.zeros(nr * nc * rescaleFac * rescaleFac * othdimprod).reshape(othdim + (nr * rescaleFac, nc * rescaleFac)) + MV
    else:
        fine = np.zeros(nr * nc * rescaleFac * rescaleFac).reshape(nr * rescaleFac, nc * rescaleFac) + MV

    dimF = np.shape(fine)
    nrF,ncF = dimF[-2:]
    ii = -1
    for i in range(0, nrF):
        if i % rescaleFac == 0:
            ii += 1

        # repeat along space axis, which (I think!!!) should
        # always be the last dimension of coarse[...,ii,]
        # hence, ndim - 2 (-1 to comply with Python zero
        # indexing, -1 because we select a single row and hence
        # reduce the number of dimensions by 1)
        fine[...,i,:] = coarse[...,ii,].repeat(rescaleFac, ndim - 2)

    # SM: not exactly sure what is going on here
    nr = None; nc = None
    del nr; del nc
    nrF = None; ncF = None
    del nrF; del ncF
    n = gc.collect(); del gc.garbage[:]; n = None; del n
    return fine

# def regridData2FinerGrid_old(rescaleFac,coarse,MV):
#     if rescaleFac ==1:
#         return coarse
#     nr,nc = np.shape(coarse)
    
#     fine= np.zeros(nr*nc*rescaleFac*rescaleFac).reshape(nr*rescaleFac,nc*rescaleFac) + MV
    
 
#     ii = -1
#     nrF,ncF = np.shape(fine)
#     for i in range(0 , nrF):
#             if i % rescaleFac == 0:
#                 ii += 1
#             fine [i,:] = coarse[ii,:].repeat(rescaleFac)

#     nr = None; nc = None
#     del nr; del nc
#     nrF = None; ncF = None
#     del nrF; del ncF
#     n = gc.collect() ; del gc.garbage[:] ; n = None ; del n
#     return fine

# def regridToCoarse(fine,fac,mode,missValue):
#     nr,nc = np.shape(fine)
#     coarse = np.zeros(nr/fac * nc / fac).reshape(nr/fac,nc/fac) + MV
#     nr,nc = np.shape(coarse)
#     for r in range(0,nr):
#         for c in range(0,nc):
#             ar = fine[r * fac : fac * (r+1),c * fac: fac * (c+1)]
#             m = np.ma.masked_values(ar,missValue)
#             if ma.count(m) == 0:
#                 coarse[r,c] = MV
#             else:
#                 if mode == 'average':
#                     coarse [r,c] = ma.average(m)
#                 elif mode == 'median': 
#                     coarse [r,c] = ma.median(m)
#                 elif mode == 'sum':
#                     coarse [r,c] = ma.sum(m)
#                 elif mode =='min':
#                     coarse [r,c] = ma.min(m)
#                 elif mode == 'max':
#                     coarse [r,c] = ma.max(m)
#     return coarse    
            
# def waterBalanceCheck(fluxesIn,fluxesOut,preStorages,endStorages,processName,PrintOnlyErrors,dateStr,threshold=1e-5,landmask=None):
#     """ Returns the water balance for a list of input, output, and storage map files  """
#     # modified by Edwin (22 Apr 2013)

#     inMap   = pcr.spatial(pcr.scalar(0.0))
#     outMap  = pcr.spatial(pcr.scalar(0.0))
#     dsMap   = pcr.spatial(pcr.scalar(0.0))
    
#     for fluxIn in fluxesIn:
#         inMap   += fluxIn
#     for fluxOut in fluxesOut:
#         outMap  += fluxOut
#     for preStorage in preStorages:
#         dsMap   += preStorage
#     for endStorage in endStorages:
#         dsMap   -= endStorage

#     a,b,c = getMinMaxMean(inMap + dsMap- outMap)
#     if abs(a) > threshold or abs(b) > threshold:
#         if PrintOnlyErrors: 
            
#             msg  = "\n"
#             msg += "\n"
#             msg  = "\n"
#             msg += "\n"
#             msg += "##############################################################################################################################################\n"
#             msg += "WARNING !!!!!!!! Water Balance Error %s Min %f Max %f Mean %f" %(processName,a,b,c)
#             msg += "\n"
#             msg += "##############################################################################################################################################\n"
#             msg += "\n"
#             msg += "\n"
#             msg += "\n"
            
#             logger.error(msg)

#             #~ pcr.report(inMap + dsMap - outMap,"wb.map")
#             #~ os.system("aguila wb.map")
            
#             #~ # for debugging:
#             #~ error = inMap + dsMap- outMap
#             #~ os.system('rm error.map')
#             #~ pcr.report(error,"error.map")
#             #~ os.system('aguila error.map')
#             #~ os.system('rm error.map')
            
#     #~ wb = inMap + dsMap - outMap
#     #~ maxWBError = pcr.cellvalue(pcr.mapmaximum(pcr.abs(wb)), 1, 1)[0]
#     #~ #return wb




# def waterBalance(  fluxesIn,  fluxesOut,  deltaStorages,  processName,   PrintOnlyErrors,  dateStr,threshold=1e-5):
#     """ Returns the water balance for a list of input, output, and storage map files and """

#     inMap = pcr.spatial(pcr.scalar(0.0))
#     dsMap = pcr.spatial(pcr.scalar(0.0))
#     outMap = pcr.spatial(pcr.scalar(0.0))
#     inflow = 0
#     outflow = 0
#     deltaS = 0
#     for fluxIn in fluxesIn:
#         inflow += getMapTotal(fluxIn)
#         inMap += fluxIn
#     for fluxOut in fluxesOut:
#         outflow += getMapTotal(fluxOut)
#         outMap += fluxOut
#     for deltaStorage in deltaStorages:
#         deltaS += getMapTotal(deltaStorage)
#         dsMap += deltaStorage

#     #if PrintOnlyErrors:
#     a,b,c = getMinMaxMean(inMap + dsMap- outMap)
#     # if abs(a) > 1e-5 or abs(b) > 1e-5:
#     # if abs(a) > 1e-4 or abs(b) > 1e-4:
#     if abs(a) > threshold or abs(b) > threshold:
#         print "WBError %s Min %f Max %f Mean %f" %(processName,a,b,c)
#     #    if abs(inflow + deltaS - outflow) > 1e-5:
#     #        print "Water balance Error for %s on %s: in = %f\tout=%f\tdeltaS=%f\tBalance=%f" \
#     #        %(processName,dateStr,inflow,outflow,deltaS,inflow + deltaS - outflow)
#     #else:
#     #   print "Water balance for %s: on %s in = %f\tout=%f\tdeltaS=%f\tBalance=%f" \
#     #        %(processName,dateStr,inflow,outflow,deltaS,inflow + deltaS - outflow)

#     wb = inMap + dsMap - outMap
#     maxWBError = pcr.cellvalue(pcr.mapmaximum(pcr.abs(wb)), 1, 1)[0]

#     #if maxWBError > 0.001 / 1000:
#         #row = 0
#         #col = 0
#         #cellID = 1
#         #troubleCell = 0

#         #print "Water balance for %s on %s: %f mm !!! " %(processName,dateStr,maxWBError * 1000)
#         #pcr.report(wb,"%s-WaterBalanceError-%s" %(processName,dateStr))

#         #npWBMError = pcr2numpy(wb, -9999)
#         #(nr, nc) = np.shape(npWBMError)
#         #for r in range(0, nr):
#             #for c in range(0, nc):

#                 ## print r,c

#                 #if npWBMError[r, c] != -9999.0:
#                     #val = npWBMError[r, c]
#                     #if math.fabs(val) > 0.0001 / 1000:

#                         ## print npWBMError[r,c]

#                         #row = r
#                         #col = c
#                         #troubleCell = cellID
#                 #cellID += 1
#         #print 'Water balance for %s on %s: %f mm row %i col %i cellID %i!!! ' % (
#             #processName,
#             #dateStr,
#             #maxWBError * 1000,
#             #row,
#             #col,
#             #troubleCell,
#             #)

#     return inMap + dsMap - outMap



# def waterAbstractionAndAllocationHighPrecision_NEEDMORETEST(water_demand_volume, \
#                                                available_water_volume, \
#                                                allocation_zones,\
#                                                zone_area = None,
#                                                debug_water_balance = True,\
#                                                extra_info_for_water_balance_reporting = ""):

#     # STILL UNDER DEVELOPMENT - NOT FULLY TESTED
    
#     logger.debug("Allocation of abstraction. - using high precision option")
    
#     # demand volume in each cell (unit: m3)
#     remainingcellVolDemand = pcr.max(0.0, water_demand_volume)
    
#     # available water volume in each cell
#     remainingCellAvlWater  = pcr.max(0.0, available_water_volume)

#     # loop from biggest values of cellAvlWater
#     min_power_number = 0                                                            # The minimum value is zero.
#     max_power_number = int(pcr.mapmaximum(pcr.log10(remainingCellAvlWater))) + 1
#     step = 1
#     cell_abstrac_for_every_power_number = {}
#     cell_allocat_for_every_power_number = {}
#     for power_number in range(max_power_number, min_power_number - step, -step):
        

#         logger.debug("Allocation of abstraction. - using high precision option - loop power number: " + str(power_number))

#         # cell available water in this loop        
#         cellAvlWater = pcr.rounddown(remainingCellAvlWater * pcr.scalar(10.**(power_number))) / pcr.scalar(10.**(power_number))
#         if power_number == min_power_number: cellAvlWater = pcr.max(0.0, remainingCellAvlWater)
        
#         # zonal available water in this loop
#         zoneAvlWater = pcr.areatotal(cellAvlWater, allocation_zones)

#         # total actual water abstraction volume in each zone/segment (unit: m3)
#         # - limited to available water
#         zoneVolDemand   = pcr.areatotal(remainingcellVolDemand, allocation_zones)
#         zoneAbstraction = pcr.min(zoneAvlWater, zoneVolDemand)
    
#         # actual water abstraction volume in each cell (unit: m3)
#         cellAbstraction = getValDivZero(\
#                           cellAvlWater, zoneAvlWater, smallNumber) * zoneAbstraction
#         cellAbstraction = pcr.min(cellAbstraction, cellAvlWater)                                                                   
        
#         # allocation water to meet water demand (unit: m3)
#         cellAllocation  = getValDivZero(\
#                           remainingcellVolDemand, zoneVolDemand, smallNumber) * zoneAbstraction 
    
#         # water balance check
#         if debug_water_balance and not isinstance(zone_area,types.NoneType):
#             waterBalanceCheck([pcr.cover(pcr.areatotal(cellAbstraction, allocation_zones)/zone_area, 0.0)],\
#                               [pcr.cover(pcr.areatotal(cellAllocation , allocation_zones)/zone_area, 0.0)],\
#                               [pcr.scalar(0.0)],\
#                               [pcr.scalar(0.0)],\
#                               'abstraction - allocation per zone/segment (with high precision) - loop (power number): ' + str(power_number) ,\
#                                True,\
#                                extra_info_for_water_balance_reporting, threshold = 1e-5)

#         # actual water abstraction and allocation in this current loop (power number)
#         cell_abstrac_for_every_power_number[str(power_number)] = cellAbstraction
#         cell_allocat_for_every_power_number[str(power_number)] = cellAllocation
                
#         # remaining cell available water and demand 
#         remainingCellAvlWater  = pcr.max(0.0, remainingCellAvlWater  - cellAbstraction)
#         remainingcellVolDemand = pcr.max(0.0, remainingcellVolDemand - cellAllocation )
        
#     # sum from the smallest values (minimizing numerical errors)
#     sumCellAbstraction = pcr.scalar(0.0)
#     sumCellAllocation  = pcr.scalar(0.0)
#     for power_number in range(min_power_number, max_power_number + step, step):
#         sumCellAbstraction += cell_abstrac_for_every_power_number[str(power_number)]
#         sumCellAllocation  += cell_allocat_for_every_power_number[str(power_number)]
    
#     # water balance check
#     if debug_water_balance and not isinstance(zone_area,types.NoneType):
#         waterBalanceCheck([pcr.cover(pcr.areatotal(sumCellAbstraction, allocation_zones)/zone_area, 0.0)],\
#                           [pcr.cover(pcr.areatotal(sumCellAllocation , allocation_zones)/zone_area, 0.0)],\
#                           [pcr.scalar(0.0)],\
#                           [pcr.scalar(0.0)],\
#                           'abstraction - allocation per zone/segment (with high precision) - sum after loop' ,\
#                            True,\
#                            extra_info_for_water_balance_reporting, threshold = 1e-5)
    
#     return sumCellAbstraction, sumCellAllocation

# def waterAbstractionAndAllocation(water_demand_volume,available_water_volume,allocation_zones,\
#                                   zone_area = None,
#                                   high_volume_treshold = 1000000.,
#                                   debug_water_balance = True,\
#                                   extra_info_for_water_balance_reporting = "",
#                                   landmask = None,
#                                   ignore_small_values = False):

#     logger.debug("Allocation of abstraction.")
    
#     # demand volume in each cell (unit: m3)
#     cellVolDemand = pcr.max(0.0, water_demand_volume)
#     if not isinstance(landmask, types.NoneType):
#         cellVolDemand = pcr.ifthen(landmask, pcr.cover(cellVolDemand, 0.0))
#     if ignore_small_values: # ignore small values to avoid runding error
#         cellVolDemand = pcr.rounddown(pcr.max(0.0, water_demand_volume))
#     else:
#         cellVolDemand = pcr.max(0.0, water_demand_volume)
    
#     # total demand volume in each zone/segment (unit: m3)
#     zoneVolDemand = pcr.areatotal(cellVolDemand, allocation_zones)
    
#     # total available water volume in each cell
#     cellAvlWater = pcr.max(0.0, available_water_volume)
#     if not isinstance(landmask, types.NoneType):
#         cellAvlWater = pcr.ifthen(landmask, pcr.cover(cellAvlWater, 0.0))
#     if ignore_small_values: # ignore small values to avoid runding error
#         cellAvlWater = pcr.rounddown(pcr.max(0.00, available_water_volume))
#     else:
#         cellAvlWater = pcr.max(0.0, available_water_volume)
    
#     # total available water volume in each zone/segment (unit: m3)
#     # - to minimize numerical errors, separating cellAvlWater 
#     if not isinstance(high_volume_treshold,types.NoneType):
#         # mask: 0 for small volumes ; 1 for large volumes (e.g. in lakes and reservoirs)
#         mask = pcr.cover(\
#                pcr.ifthen(cellAvlWater > high_volume_treshold, pcr.boolean(1)), pcr.boolean(0))
#         zoneAvlWater  = pcr.areatotal(
#                         pcr.ifthenelse(mask, 0.0, cellAvlWater), allocation_zones)
#         zoneAvlWater += pcr.areatotal(                
#                         pcr.ifthenelse(mask, cellAvlWater, 0.0), allocation_zones)
#     else:
#         zoneAvlWater  = pcr.areatotal(cellAvlWater, allocation_zones)
    
#     # total actual water abstraction volume in each zone/segment (unit: m3)
#     # - limited to available water
#     zoneAbstraction = pcr.min(zoneAvlWater, zoneVolDemand)
    
#     # actual water abstraction volume in each cell (unit: m3)
#     cellAbstraction = getValDivZero(\
#                       cellAvlWater, zoneAvlWater, smallNumber)*zoneAbstraction
#     cellAbstraction = pcr.min(cellAbstraction, cellAvlWater)                                                                   
#     if ignore_small_values: # ignore small values to avoid runding error
#         cellAbstraction = pcr.rounddown(pcr.max(0.00, cellAbstraction))
#     # to minimize numerical errors, separating cellAbstraction 
#     if not isinstance(high_volume_treshold,types.NoneType):
#         # mask: 0 for small volumes ; 1 for large volumes (e.g. in lakes and reservoirs)
#         mask = pcr.cover(\
#                pcr.ifthen(cellAbstraction > high_volume_treshold, pcr.boolean(1)), pcr.boolean(0))
#         zoneAbstraction  = pcr.areatotal(
#                            pcr.ifthenelse(mask, 0.0, cellAbstraction), allocation_zones)
#         zoneAbstraction += pcr.areatotal(                
#                            pcr.ifthenelse(mask, cellAbstraction, 0.0), allocation_zones)
#     else:
#         zoneAbstraction  = pcr.areatotal(cellAbstraction, allocation_zones)    
    
#     # allocation water to meet water demand (unit: m3)
#     cellAllocation  = getValDivZero(\
#                       cellVolDemand, zoneVolDemand, smallNumber)*zoneAbstraction 
    
#     # extraAbstraction to minimize numerical errors:
#     zoneDeficitAbstraction = pcr.max(0.0,\
#                                      pcr.areatotal(cellAllocation , allocation_zones) -\
#                                      pcr.areatotal(cellAbstraction, allocation_zones))
#     remainingCellAvlWater = pcr.max(0.0, cellAvlWater - cellAbstraction)
#     cellAbstraction      += zoneDeficitAbstraction * getValDivZero(\
#                             remainingCellAvlWater, 
#                             pcr.areatotal(remainingCellAvlWater, allocation_zones), 
#                             smallNumber)                        
#     # 
#     # extraAllocation to minimize numerical errors:
#     zoneDeficitAllocation = pcr.max(0.0,\
#                                     pcr.areatotal(cellAbstraction, allocation_zones) -\
#                                     pcr.areatotal(cellAllocation , allocation_zones))
#     remainingCellDemand = pcr.max(0.0, cellVolDemand - cellAllocation)
#     cellAllocation     += zoneDeficitAllocation * getValDivZero(\
#                           remainingCellDemand, 
#                           pcr.areatotal(remainingCellDemand, allocation_zones), 
#                           smallNumber)                        
    
#     if debug_water_balance and not isinstance(zone_area,types.NoneType):

#         waterBalanceCheck([pcr.cover(pcr.areatotal(cellAbstraction, allocation_zones)/zone_area, 0.0)],\
#                           [pcr.cover(pcr.areatotal(cellAllocation , allocation_zones)/zone_area, 0.0)],\
#                           [pcr.scalar(0.0)],\
#                           [pcr.scalar(0.0)],\
#                           'abstraction - allocation per zone/segment (PS: Error here may be caused by rounding error.)' ,\
#                            True,\
#                            extra_info_for_water_balance_reporting,threshold=1e-4)
    
#     return cellAbstraction, cellAllocation

def findLastYearInNCFile(ncFile):

    # open a netcdf file:
    if ncFile in filecache.keys():
        f = filecache[ncFile]
    else:
        f = nc.Dataset(ncFile)
        filecache[ncFile] = f

    # last datetime
    last_datetime_year = findLastYearInNCTime(f.variables['time']) 
    
    return last_datetime_year
    
def findLastYearInNCTime(ncTimeVariable):

    # last datetime
    last_datetime = nc.num2date(ncTimeVariable[len(ncTimeVariable) - 1],\
                                ncTimeVariable.units,\
                                ncTimeVariable.calendar) 
    
    return last_datetime.year

def findFirstYearInNCTime(ncTimeVariable):

    # first datetime
    first_datetime = nc.num2date(ncTimeVariable[0],\
                                ncTimeVariable.units,\
                                ncTimeVariable.calendar) 
    
    return first_datetime.year

def cmd_line(command_line,using_subprocess = True):

    msg = "Call: "+str(command_line)
    logger.debug(msg)
    
    co = command_line
    if using_subprocess:
        cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open('/dev/null'),shell=True).communicate()
    else:
        os.system(co)

# def plot_variable(pcr_variable, filename = "test.map"):

#     pcr.report(pcr_variable, filename)
#     cmd = 'aguila '+str(filename)
#     os.system(cmd)
    
    
